
#https://www.youtube.com/watch?v=QljEeBei-JA Vegan R Library Tutorial

library(dplyr)
library(vegan)
library(ggplot2)
library(ggthemes)
library(ggrepel)
library(shadowtext)
library(ggforce)


datos1 <- read.delim('datos_limpios.csv', sep = ';', header = TRUE)


# Age
datos1$cl_age_interval <- factor(datos1$cl_age_interval,
                                 labels = c('15-24','25-34','35-44','45-54',
                                            '55-64','65+'))
####
## Sample Descriptive Statistics Table
####

age_tbl <- as.data.frame.matrix(table(datos1$cl_age_interval, datos1$cl_sex_factor))

## Attitude toward AI, age and sex
datos <- read.delim('age_sex_t_AI.csv', sep = ';', header = TRUE, fileEncoding="latin1")

#colnames(datos)[4:8] <- c('Totally \ndisagree', 'Tend to \ndisagree', 'Neither agree \nnor disagree',
#                          'Tend to \nagree','Totally \nagree')
colnames(datos)[4:8] <- c('VN', 'N', 'NE','P','VP')
row.names(datos) <- datos$cl_group

aa <- datos1 %>%
  group_by(cl_age_interval, cl_sex_factor) %>%
  summarise(across(everything(),mean),
            .groups = 'drop')

a <- aa %>% 
  filter(., cl_sex_factor == 'Woman')
b <- aa %>% 
  filter(., cl_sex_factor == 'Man')

datos.means <- bind_rows(a,b)

###### KEY TO SELECT THE CORRECT VARIABLES
model <- select(datos, -starts_with('cl'))
# model1 <- select(datos.means, starts_with(c('i','g')))
model1 <- select(datos.means, starts_with(c('in_')))

#ordination by NMDS
NMDS <- metaMDS(model, distance = "bray", k = 2)
## ggplot
nmds.plot <- data.frame(name = rownames(model), axis1 = NMDS$points[,1], axis2 = NMDS$points[,2])
nmds.species <- data.frame(name = colnames(model), axis1 = NMDS$species[,1], axis2 = NMDS$species[,2])

stressplot(NMDS)

plot(NMDS)

##### Environmental variables
ef <- envfit(NMDS, model1[2:ncol(model1)])
ef

#Get the vectors for env.fit
df_envfit<-scores(ef,display=c("vectors"))
df_envfit<-df_envfit*vegan:::ordiArrowMul(df_envfit)
#size adjusted by R2
df_envfit<-as.data.frame(df_envfit)*ef$vectors$r
#
#add p-values
df_envfit <- cbind(df_envfit, data.frame(pvals = ef$vectors$pvals))


#Only significant, less th
######
df_envfit <- dplyr::filter(df_envfit, pvals < 0.1)

#####
gof <- goodness(NMDS)
gof

nmds.plot$grp<-rep(c("Woman","Man"),each=6)
nmds.plot$gof <- 2*gof/mean(gof)
nmds.species$color <- c('n', 'n', 'ne', 'p', 'p')

## NEW Change Format Plot
#png(file="plot1.png", width=400, height=400)

ggplot() +
  geom_text_repel(data=as.data.frame(df_envfit*0.6),
                  mapping = aes(NMDS1 * 1.95, NMDS2 * 1.95, label = rownames(df_envfit)),  size=4, ## Environment
                  color="black", alpha = 0.9, inherit.aes = FALSE) +
  geom_text(nmds.plot, mapping = aes(axis1, axis2, label = name, color = grp),
            alpha = 0.9, inherit.aes = FALSE,  fontface = "bold", size = 4) + ## Sites
  theme(legend.title=element_blank(),legend.position = 'none') +
  geom_point(nmds.plot, mapping = aes(axis1, axis2, label = name, color  = grp, size=1*gof), ## Circle
             inherit.aes = FALSE, shape = 21, fill = NA, stroke = 1, alpha = 0.6) +
  scale_size_continuous(range = c(3, 12)) +
  geom_text(data = nmds.species, mapping = aes(axis1, axis2, label = name, color = color),
            color = c("firebrick3", "firebrick3", "gray48","springgreen3","springgreen3"),
            inherit.aes = FALSE,  alpha = .6,  size=6, fontface = "bold") + ## Species
  geom_segment(data=df_envfit, mapping = aes(x = 0, y = 0, xend = NMDS1*0.8, yend = NMDS2*0.8), ## Arrow
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")),color="black",
  alpha=0.9, lineend = 'round', linejoin = 'round', size = 0.7) +
  # xlim(-0.20, 0.20) +
  # ylim(-0.10, 0.10) +
  geom_hline(yintercept = 0, alpha = 0.2) +
  geom_vline(xintercept = 0, alpha = 0.2) +
  labs(x = '', y = '', color = NULL) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), panel.border = element_blank()) +
  guides(x = "none", y = "none")



#####################
#Bootstrapping and testing for differences between the groups, should be small, the groups are very different from each other
fit <- adonis2(model ~ cl_age + cl_sex , data=datos, permutations=9999, method="bray")
fit


#####################
#Check assumption of homogeneity of multivariate dispersion P should be > 0.05 to have no problems with the assumption
distances_data <- vegdist(model)
anova(betadisper(distances_data, datos$cl_sex))
anova(betadisper(distances_data, datos$cl_age))
anova(betadisper(distances_data, datos$cl_group))
### Stress
NMDS[["stress"]]

